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ABSTRACT 



A three dimensional staggered grid model, based on the Boussinesq 
equations, is developed, and applied to mesoscale waves. The model is 
tested with a three dimensional stable gravity wave and forecasts for 
40 time steps are compared with the analytical solutions and also with 
a full grid model. Comparisons show very good agreement with both 
forecasts for the stated time. Failure of the relaxation scheme to 
meet the convergence criterion after 47 time steps is discussed. 
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I. INTRODUCTION 



In the past few years many attempts have been made to understand 
the complex dynamic interactions between the thunderstorm and its 
environment. Newton and Newton (1959) proposed a dynamical thunder- 
storm model which proposed that strong vertical motions within the 
thunderstorm cell create an effective barrier to the flow of the 
environmental air. 

Ogura (1963) first discussed the probability of using the primitive 
equations for forecasting small scale dynamical features. Interaction 
between the convective elements and the environmental wind field can 
be investigated numerically by the use of these equations. Deardorff 
(1970) made a numerical study of turbulent flow using a three-dimensional 
primitive equation model. A similar forecasting procedure is used in 
this model. The importance of gravity waves in small scale convective 
activity is discussed by Matsumoto and Ninomiya (1969). Orszag (1969) 
presented a model for the simulation of turbulence which employs a 
fast Fourier Transform instead of the Liebmann Relaxation Technique. 
Further refinement of this model may be accomplished by the use of 
such a method. 

Norton (1970) discussed the development of a three dimensional model 
for application to mesoscale motions based on the Boussinesq equations. 

The model presented here is a further development of the same model with 
a staggered grid system applied to minimize computer storage requirements. 
The total grid size is 14 by 14 in the horizontal with 10 levels in the 
vertical, while the computational grid size is 11 by 10 by 6. The grid 
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interval in all directions is 1 kilometer. The Coriolis parameter, 
treated as a constant, is determined at a mean latitude of 35° North. 

Numerical computations are made using the IBM 0S/360 Computer of 
the W. R. Church Computer Center. To conform to computational stability 
a time step of 30 seconds was chosen. 



II. MODEL DEVELOPMENT EQUATIONS 



A. FORECAST EQUATIONS 

The pressure term in the complete equation of motion can be written 



as: 



- v, p 
o 3 



ir ’3 

P 3 



,1/k 

rp 



>1/k 



or 



p v 3 P = 0 V 3 8 



where 



6 = C r 



£_ 



1 K 



and 



0 = T 



P J 



Using these results in the component form of the equation of motion 
yields : 



du 

dt 


- 0 


91 + 
9X 


fv 


+ vV^ U , 


(1) 


Q-lQ. 

c+|< 

II 


- 0 


36 

3y 


fu 


2 

+ vV^ v , 


(2) 


dw 

dt 


- 0 


36 

3y 


g 


2 

+ VVg W 


(3) 



The perturbation expressions can be formed by: 
6 = B 0 (z) + 6' (x,y,z,t) 

e = 0 q + e(x,y,z,t) 

where 



|6‘ | « 6 and |e| « q q 

and the latter inequality corresponds to the Boussinesq approximation. 
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After substitution into equations (1), (2), and (3), the equations 



of motion can be written: 





du 

dt 


■ 0 o 


3(3 1 

3X 


+ fv 


+ 


vV 3 


U, 


(4) 




dv 

dt 


- 0 o 


33 1 

ay 


- fu 


+ 


vV 3 


V, 


(5) 


and 


dw 

dt 


- 0 o 


ffo 

0Z 


- 6 o 


36* 

9Z 


- -0 


36 o + > 

— - g +w 3 w . 





with the relatively small products: 



a 3JL fl Ml a nd a liLl 

6 ax 5 0 ay ’ and 0 az 



being neglected. The hydrostatic equation relates 3 0 and e Q by: 



ffc 

3Z 



S_ 



and the vertical component can now be written: 



dw 

dt 



9 |il + ae +VV 2 W 

0 9Z 0^ 3 



( 6 ) 



By defining <j> = e o e 1 and substituting into (4), (5), and (6) and 
expanding the total derivatives, the system of equations becomes: 



where 



3U 

3t 


■ - L < u > - & * * 


2 

+ vV 3 U, 


(7) 


3 V 
3t 


* -L(v) - § - fu 


2 

+ vV 3 V, 
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3W 

3t 
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+ vV 3 W, 


(9) 


39 

3t 


= -L(e) 
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(10) 


V 3 * 


V = 0 
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(ID 


L(S) 


= -v 3 -(v 3 s). 


S = U , V ,W , 0 . 
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Incompressible flow is assumed, to eliminate sound waves. This system 
has become known as the anelastic equations and are quite accurate for 
shallow convection. (Ogura and Phillips, 1962) 



B. DETERMINATION OF THE PRESSURE TERMS 

Equations (7), (8), and (9) may be combined in vector form and 
wri tten: 

3V", _ nfi _ ~ 

— ' - (V 3 ' V)V 3 - v 3* - f(k x V 3> + ^ k + v v 3 V 3 <’2) 



By taking the three dimensional divergence and assuming the time rate of 
change to be zero, (12) becomes: 




- -'3 [( V 3-'3 )V 3 ] - V 3* 

-v 3 -[f(F x V 3 )] +f- ff + Y v7 3 ¥ 



03) 



For the horizontal extent used in this model, the Coriolis parameter may 
be treated as a constant and then (13) becomes: 

- -v 3 .[(V 3 .v 3 )V 3 ] - f 0 v 3 -(k X v 3 ) + f- H 

= -» 3 -[(V 3 -v 3 )V 3 ] + f 0 ? + f- H +v 3 .w 2 V 04) 

0 

Note that when the terms on the RHS are solved by ordinary differencing 
techniques, the value of $ on the LHS can be solved by using the 
Extrapolated Liebmann Relaxation Scheme. 
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III. FINITE DIFFERENCING 



Numerical integration of equations (7) through (10) is carried out 
using finite difference methods in x,y, z, and t. Extrapolated Liebmann 
Relaxation is used to solve for the pressure term in equation (14). For 
the complete development of the relaxation scheme, see Appendix A. The 
technique applied to generate the forcing function used in (14) is shown 
in complete form in Appendix B. 

To maintain consistency in finite differencing throughout the fore- 
cast procedure, the differencing scheme used in equation (14) must be 
equivalent to the type used on equations (7), (8), and (9). The flux 
form differencing method devised by Arakawa (1966) was utilized in the 
form: 



L(s) ij, k ' 



( U i+l,j,k + u i,j,k^ ^ S i+l,j,k + S i,j,k^ 



-(u i-l,j,k + u i,j,k } (S i-l,j,k + S i,j,k ) 



^ V i >j+l >k +V i,j,k^ ( S i,j+l,k + S i,j,k^ 



/ 4ax 



~^ V i,j-l,k + v i,j,k^ ^ S i,j-l,k + S i,j,k^ 



^ w i,j,k+l W i , j ,k^ ^ S i,j,k+1 +S i,j,k^ 



/ 4&y 



~^ W i»j.k-l + w i,j,k^ ^ S i,j,k-1 + S i,j,k^ 



/ 4az 



This method is used for the u,v,w and 0 fields, but for the linear space 
derivatives for <j>, an averaging technique must be employed due to the 
staggered grid system. The centered form of the space derivative takes 
on the form: 
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(S. . . 4 S. . . , 4 s. . , . 4 s. . , , , ) / 4. 

v i,J>k i,j,k-l i,j-l,k i,j-l,k-l ; ' 

~ (S i-l,j,k + S i-l,j,k-l + S i-1 ,j-l ,k + S i-l,j-l,k-l 




/AX 



with similar forms for the y and z derivatives. The primed subscripts 
refer to the grid points one-half grid distance above or below the 
diagonals connecting surrounding velocity points in a horizontal plane. 
Centered time differencing ("leapfrog") is used for all forecasts with 
the exception of the first time step where a single forward time step 
is used. 
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IV. BOUNDARY CONDITIONS 



The boundary conditions are selected to enforce mass and energy 
conservation within the confines of the computational grid. Figure 1 
illustrates the horizontal grid domain and boundary placement. Cyclic 
continuity of parameters is imposed by forcing periodicity in the 
east-west direction according to the relations 

5 1 , J,K = S IM-2,J ,K 

5 2, J ,K = S IM- 1 ,J ,K 
S IM,J,K = S 3,J,K 

where S = u,v,w, e. Since values of S are required at three consecutive 
points in solving for the pressure and divergence terms, values outside 
the domain can be obtained through the cyclic conditions. The east- 
west boundary conditions on <j> are also cyclic and are specified by: 

<{> 1,J,K = ^IM-2 ,J ,K 

♦iM-l.J.K = ^2,J ,K . 

The south boundary is located halfway between velocity grid points 
J = 2 and J = 3, while the north boundary is between the points 
J = JM - 2 and J = JM - 1 . To prevent the flow of mass and energy 
through these boundaries, the v-component is set by the following 
conditions : 



1,1, K = 


" V I,4,K 


1,2, K = 


~ V I ,3 ,K 
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. Illustration of the grid system in the horizontal 
direction. 
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V I,JM,K " V I,JM-3,K 

V I , JM-1 , K = _V I,JM-2,K . 
The values of u,w, and e are extrapolated outward by: 



S I,1,K " S I,2,K = S I ,3,K 
S I ,JM,K = S I,JM-1,K = S I , JM-2 ,K 



where S = u,w,9. This implies that there is no advection through this 
boundary and furthermore no heat exchange as discussed by Ukaji and 
Matsuno (1970). 

The geostrophic relationship: 




was used to obtain the values of <j> outside the north-south boundaries. 

★ 

Where u • , m represents the volume average of the u velocity field at 

I 9 J 9& 

the boundary. The results reported here were for the special case of 
f = 0. For this case, the boundary condition then becomes: 

I = 0 . 

i 1 » j 1 X 




Figure 2 illustrates the vertical grid domain and boundary placements. 
The top and bottom boundaries are set in manner similar to the north- 
south boundaries where the points midway between K = KM - 1 and K = KM - 2 
and K = 2 and K = 3 are the top and bottom boundaries respectively. The 
no-flux condition is imposed on the w field by: 



I,J,1 = 




I >J >2 = 
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Figure 2. Illustration of the grid system in the vertical 
direction. 
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W I,J,KM 






W I,J,KM-1 


" W I,J,KM-2 , 


while the 


e values are set by the 


relations: 




9 I,J,1 






0 I,J,2 






0 I ,J ,KM 


-9 I ,J ,KM-3 




6 I ,J ,KM-1 


" 0 I ,J ,KM-2 . 


The values 


of u and v are simply extrapolated outward by: 




S I,J,1 


S I , J ,2 S I , J ,3 




S I ,J ,KM 


S I ,J ,KM-1 = S I,J,KM-2 



where S = u,v. 

The value of the pressure term is obtained by the use of: 




The importance of proper boundary conditions cannot be over- 
emphasized. This model was found to be extremely sensitive to all the 
boundary conditions. The complete development of the boundary conditions 
is found in Appendix C. 
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V. INITIAL CONDITIONS 



Initialization of the u,v,w and 0 fields is accomplished by the use 
of linearized gravity wave solutions, which are exact for small amplitude 
gravity waves and are derived from the following equations: 

3u _ 3 4 > 

3t " 3x 



3V_ _ 3j> 

3t “ " 3y 

3W _ 3 <f> 4 q9 

3t “ “ 3Z 0 

0 

v 3 • V = 0 
and 

li = _ w 

3t 3z 



The exact solutions for a stable wave propagating in the x direction are 
given by: 

AmW / \ 

u = — 2 2 sm(A x -ct) cos yy cos mz , 

A + y 

v _ — - cos (^ x - c t) sin yy cos mz , 

A + y 

w = W cos (a x -ct) cos yy sin mz , 

W 90 

6 = Ac 37 x -ct ) cos ^y s '’ n mz • 

where c represents the phase speed of the stable gravity wave in m sec - ^ 
and is calculated from: 



c 



— J * 2 * 

X \/,2 . 2 . J. 

"A + y + m 



( 15 ) 
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and N is the Brunt Vaisala frequency 



N 



= Js— M. 

Vo dz 



The initial guess for the 4 field is chosen to be an arbitrary constant 
(2.0) due to the extremely efficient relaxation technique. The set of 
exact solutions with t = 0 was used to initialize the fields. 
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VI . RESULTS 



By introducing a stable gravity wave perturbation into the basic 
model, the wave form can be tested for proper phase speed, amplitude 
growth or decay and computational stability. The initial conditions 
should result in a wave that propagates in the x direction alone with 
no growth in amplitude. Figures 3a, 3b, 3c, and 3d are representative 
plotted values of the u,v,w, and e wave forms after time steps 10, 

20, 30, and 40, with the initial value present for comparison. These 
plots are for the wave forms in the x direction for J = 5 and K = 7 
and show the phase speed relationship between the various time steps. 

The exact solution to equation (15) yields a phase speed of 5.5 m sec"^ 

-1 -3 -1 

for the specified values W = .01 m sec and N = 4.04 x 10 sec . 
Phase speed calculated from the plotted values is 5.4 m sec -1 for each 
of the four wave forms. There is negligible amplitude growth or decay 
in any of the plots. 

Figures 4a, 4b, 4c, and 4d are the plotted values of the wave forms 
in the y direction for 1=5 and K = 7. The plots indicate there is 
no propagation in the y direction. Due to the nature of the solutions 
to the linearized equations, the wave forms should oscillate about the 
zero value as indicated in the stated figures. 

Plots of the wave forms in the z direction for 1=5 and J = 5 are 
shown in Figures 5a, 5b, 5c, and 5d and are similar in amplitude 
oscillation to the y direction with no phase speed observed. 

The main objective of this staggered grid model was to reduce the 
storage size and computer time required for an operational small scale 
convective model. The previous full-grid model required 502K bytes of 
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storage and approximately 20 minutes of computer time to produce a 40 
time step forecast. The forecasts produced by the staggered grid 
model accomplishes a 40 time step prognosis in less than 10 minutes of 
computer time and requires only 156K bytes of storage. 

The results obtained after 40 time steps compared to the analytical 
solutions illustrate the accuracy of the finite difference solutions. 
However, a few time steps later, the solutions fail to converge. A 
further discussion of this difficulty will be given in a later section. 
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Figure 3a. u component of stable gravity wave in x direction for J 
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Figure 3b. v component of stable gravity wave 
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Figure 3c. w component of stable gravity wave in x direction for 
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Figure 3d. e wave in x direction for J 
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Figure 4a. u component of stable gravity wave in y direction for 
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Figure 4b. v component of stable gravity wave in y direction for 
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Figure 4c. w component of stable gravity wave in y direction for 
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Figure 4d. e wave in the y direction for 
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Figure 5a. u component of stable gravity wave in z direction for 
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Figure 5b. v component of stable gravity wave in z direction for 
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Figure 5d. e wave in z direction for 



VII. CONCLUSIONS AND SUMMARY 



A ten-level staggered grid primitive equation model was applied to a 
mesoscale grid and tested with an initial stable gravity wave perturba- 
tion to check the computational stability and finite difference form of 
the equations of motion. Comparing the results presented here to a non- 
staggered grid model, the computer time and storage requirements have 
been reduced to about 1/3 of the full grid model. 

The forecasts produced by this model compare quite closely with the 
full grid model but ceases to function after 47 time steps. Up until 
this time step all fields appear to be predicted properly. Tests were 
run in double precision to determine if truncation error was a problem 
and the same type of results occured with the model failing to converge 
after 47 time steps. In all cases the points that failed to converge 
were near the boundaries which leads to the conclusion that further 
refinements of the boundary conditions may be necessary. Experiments 
have shown that a less stringent relaxation tolerance produces less 
exact <f> gradients which in turn induces greater errors in the prognosis 
of the other parameters. Apparently there is still some difficulty in 
the program and presently studies are being co. ducted to further refine 
this model . 

Continued studies using this model should have a provision for 
including vertical compressibility. Additionally, since a primary 
source of the energy in a thunderstorm- type cloud is latent heat, 
application to cumulus scale convective activity would be inadequate 
without the inclusion of moisture. 



36 



APPENDIX A 



EXTRAPOLATED LIEBMANN RELAXATION SCHEME 



FOR THREE DIMENSIONAL APPLICATION 



The equation to be solved by relaxation is: 




(A-l ) 



where the tendency of the divergence is zero and: 



F = -v 3 - (V 3 • v 3 )V 3 + f 0 5 + 



3_ ii 

e o 3Z 



(A- 2 ) 



with the <}> points staggered in the vertical as well as in the horizontal. 
(They are located one-half grid distance above or below the diagonals 
connecting velocity points as in Fig. Al.) In other words, the velocity 
and e points are located at the corners of a cube and the <f> point in the 
center. For convenience, a new temporary indexing scheme will be 
util ized by writing: 



where i,j,k represent u,v,w and 0 points and i 1 , j 1 , k' represent <|> 
points. Writing (A-l) with this new notation in finite difference form 
at the point i ' , j', k' and expanding yields: 



• I 



2 



j 1 = j - 2 



k 



k 



2 
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Figure A1 . Three dimensional schematic of the location of 
Various fields. 
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(4a y)‘ 



+4 V-1 


J'-l 


,k'+l 


+ 2 *1'+l,j'+l,k' + 


4 i' 


+lj' 


,k' + 


2 V+1 


- 4 *,, 


j+1 ,k 


i -84) 


i ' , j ' ,k' ' 4<J) i 1 ,j-l 


> k ' 


+ 24 


i '-1 ,j 


' +1 , k ' 


+4 4 > i i _ 


l.j'. 


k' + 


2<f) i '-1 ,j '-1 , k ' + 


'+1, 


j'+l 


»k'-l 


+ 2<|> . , 


+<f> i'+l 


.j'-l 


»k'-l 


" 24> i * ,j '+1 ,k'-l " 


44. 


' »j' 


»k'-l 


- 2+,.. 


+4> i *-l 


.j'+l 


*k'-l 


+ 2< V-l,j',k'-l + 




-1 ,j 


'-1 , k 1 


-]* 


*i'+l. 


j'+l » 


k'+l 


" 2 *i'+l,j,k'+l + * 


i'+l 


» j ' - 


1 , k'+l 


+ 2<f>.j 



(4az)‘ 



" 44> i' ,j' ,k'+l + 2 *i' , j'-l, k'+l + V-l, j'+l, k'+l ' 2 *1 '-1 ,j ' ,k'+l 
+ * i'-l, j'-l, k'+l + 2 *i'+l,j'+l,k' " 4< N '+1 ,j 1 ,k' + 2<t> i '+1 ,j '-1 ,k' 
+4<J) i ' , j '+1 , k ' ' 8<}, i',j',k' + 4 *i',j'-l,k' + 2 *i'-l,j'+l,k' 

” 4< N ' - 1 , j 1 , k ' + 2< h '-1 ,j '-1 ,k' + '+l,j'+l,k'-l ' 2 *i'+l,j',k'-l 

+ V+l,j'-l,k'-l + 2<f> i ' ,j '+1 ,k'-l " 4< V,j',k'-l + 2 *1 1 ,j'-l ,k'-l 

+4> i ’ - 1 , j '+1 ,k'-l ' 2(f, i'-l,j' ,k'-l + *i'-l,j'-l,k'-l 

*i'+l ,j'+l ,k'+l + 2< N '+1 ,j 1 k'+l + V+l,j'-l,k'+l + 2<j) i 1 ,j '+1 ,k'+l 

+ 4 <f, i',j',k'+l + 2 *i 1 , j'-l ,k'+l + *i'-l, j'+l, k'+l + 2< N '-1 ,j ' ,k'+l 
+< ^i '-1 ,j'-l ,k'+l " 2< N '+1 , j '+1 ,k' " 4 < N'+l,j',k' " 2< * > i '+1 ,j '-1 ,k' 
" 4 *i',j'+l,k' ‘ 8 *i 1 ,j* ,k' ‘ 4 < N',j'-l,k' " '-1 ,j '+1 ,k' 

" 4<t> i ' -1 , j ' , k ' " 2 *i'-l,j'-l,k' + V+l,j'+l,k'-l + 2< N '+1 ,j ’ , k ' - 1 
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+ 4 > 



i'+l.j'-l.k'-l + 2 *i',j'+l,k'-l + 4<f, i' ,j' ,k'-l + 2<f, i' .j'-l.k'-l 



^i'-lj'+l.k'-l + 2 *i'-l,j' ,k'-l + *i'-l, j'-l.k'-l 



" F i ' »j 1 »k' = 0 



(A-3) 



Rearranging terms and clearing fractions leads to: 



(4Ay) 2 (4 az) 2 



+(4ax) 2 (4az) 2 



+(4ax) 2 ( 4Ay ) 2 



B " 8 *i',j'.k' 



C ‘ 8 *i',j',k'_ 
D " 8<f) i ' ,j ' ,k' 



-(4ax) 2 ( 4Ay ) 2 (4az) 2 F., = 0 

• * J > N 



(A-4) 



where B, C and D represent the finite difference forms in equation (A-3) 

with the term 8(f)., , subtracted off. With further rearranging of 

1 5 J > K 

terms ^\-4) becomes: 



♦i'.j'.k' = A 



(4Ay) 2 (4 az) 2 B + (4ax) 2 (4az) 2 C 



+ (4ax) 2 ( 4Ay ) 2 D - (4ax) 2 (4Ay) 2 (4 az) 2 F. , , 

• J J > ^ 



where: 



A = 



8 |( 4Ay ) 2 (4 az) 2 + (4ax) 2 (4az) 2 + (4ax) 2 (4Ay) 2 
old 



Adding and subtracting <}>., , yields: 

1 > J > K 



new _ A old .. 
*i'.j'.k' " *i ' , j ' ,k 1+aA 



(4Ay) 2 (4 az) 2 B + (4ax) 2 (4az) 2 C + 



+(4ax) 2 ( 4Ay ) 2 D - (4ax) 2 (4Ay) 2 (4 az) 2 F-, ,, 

1 yJ y^ 



(A- 5) 
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Where the over-relaxation coefficient a is introduced. (A-5)then 
becomes : 



*1\j',k' = (1 " a) ♦i'J'.k' + a A 



old 



(4Ay) 2 (4 az) 2 B 



+ (4ax) 2 (4az) 2 C + (4ax) 2 (4Ay) 2 D 



- (4ax) 2 (4Ay) 2 (4 az) 2 F., J1>k , 



(A-6) 



This was the form used for relaxation in the model. Many tests were run 
to determine the most efficient relaxation coefficient for this type of 
finite difference scheme. The value of a which resulted in the most 
efficient relaxation was 1.05. Note that if ax = Ay = az, equation (A-6) 
reduces to the more familiar form: 



.new _ n » .old 

*i\j',k' " (1 ‘ a) V 



j ' ,k' 24 



B + C + D 



‘ (4Az)2 F i ' , j 1 ,k' 

2 

The stencil of the method of the | operator in the three dimension- 
al staggered grid, assuming that the three grid lengths are equal, is 
shown in Fig. A2. 
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Figure A2. f 2 Operator in three dimensional form assuming 

AX = Ay = AZ. 
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APPENDIX B 



DEVELOPMENT OF THE FORCING FUNCTION 



The formation of the forcing function used in the relaxation scheme 
is a long process. Consistency between the finite difference methods 
used in the prognostic equations and in the relaxation technique must 
be maintained. The computation of the forcing function given in 
equation (4-2) in finite difference form follows. In order to determine 
this forcing function in a consistent manner, the forcing function must 
first be calculated at an i,j,k point and then volume averaged to 
have the function at an i',j',k' point. Each of the three directional 
components are computed separately and then summed to put the complete 
forcing function at the desired point. For example, in order to form 
F 5 , 5 , (the forcing function at the point i 1 = 5, j 1 =5, k' =5) 
the following steps are required: 





I 



1 _ 

AX 







AZ 



F 




+ F Z 

h 5',5',5 
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V 
















where 



and the 



G(u) = - L(u) + fv + vVg u 

G(v) = - L(v) - fu + vV^ v 

G(w) = - L(w) + + vV? w 

9 o d 

functions are calculated on the i, 



i,k grid. 
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APPENDIX C 



COMPLETE DEVELOPMENT OF THE BOUNDARY CONDITIONS 



For the boundary described by the plane of j ' 
of the equation of motion may be written 



2* the y component 



9V 

9t i',2,k' 



= 0 



which implies 



3 v 

3t. 



3 V 
3t • 



'i,3,k 0 °i,2,k 

Further expansion of each of the terms leads to: 



(C-l) 



and 



!¥. = _ I ( v ). - M. - f u. 

3 1 i ^ i j 3 , k 3y. ^ o i » 3 5 k 



"ft “ L ( v )i 2 k + fy" + f o u i 2 k 
9t i ,2,k 1j ^’ K 9y i ,2 , k 0 1 »^» lc 



(C-2 



(C-3) 



By introducing the no-flux conditions on v at the lateral walls from 
(C-l), (C-2) and (C-3) yield: 



_ -fu = M. +fu 

ay i,3,k 0 i,3,k 9y i,2,k 0 1,2,k 

Expansion of (C-4) into finite differences and averaging the u 
field to be consistent with the field yields: 

♦r,l,k' ^i 1 ,3,k' +A 2T 1 (u i,2,k + u i+l,2,k + u i,2,k+l 



+ u i+l ,2,k+l ^ 



(C-4) 
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with the no-flux conditions being: 



v i,2,k = " v i ,3,k 



' " v i ,4,k 

A similar approach is taken at the northern boundary with the following 
conditions : 

v i ,J1 ,k = ~ v i ,J2,k 



v i ,JM,k " " v i,J3,k 



V.JT.k' = <f> i\J3,k' " ^ u i ,J2,k + u i+l,J2,k 



+u i , J 2 , k+1 + u i+1,J2,k+l ^ 



The case considered here was for f = 0. 

For the boundary described by k' =2', the z component of the 
equation of motion may be written: 



aw 

3t i 1 ,2' , k ' 



= 0 



which implies: 



aw 

at i,3,k 



aw 

3t i,2,k 



(C-5) 



Further expansion of each of the terms in (C-5) yields: 



fr - - i_(w). 3 k - H 

3t i ,3,k 1,J ’ K 3 i,3,k e o 



aw 

at 



0 . 



i ,3,k 



l.2.k'* L<W> i' 2 -k + »,. 2 ,k'^ 6i ' 2 ' k 



(C-6) 

(C-7) 
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By introducing the no- flux conditions on w at the top and bottom 



boundaries from (C-5), (C-6) and (C-7) they yield: 



9<j> 

az 



0 . 



i >3,k 



= . a_ 

i.B, k az i>2jk 0 O 



0. 



i,2,k 



(C-8) 



By setting 0 at the boundaries as previously discussed, (C-8) is 
satisfied if: 



V 



,jM' 



= <P 



i 1 > j 1 >3 ' 



where the 0 boundary conditions are: 



0 i»j,2 = " 0 i,j,3 
0 i >j»l " " 0 i.j>4 . 



A similar approach is taken at the top boundary with the following 
conditions : 

0 i,j,Kl = -9 i , j ,K2 

9 i > j >KM = ~ 0 i , j ,K3 

and 

9 i 1 s j 1 »K1 1 = 9 i 1 , j 1 , K3 1 . 
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